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ABSTRACT 

Context. In previous work, we developed a quasi-Gaussian approximation for the likelihood of correlation functions that incorporates 
fundamental mathematical constraints on correlation functions, in contrast to the usual Gaussian approach. The analytical computation 
of these constraints is only feasible in the case of correlation functions of one-dimensional random fields. 

Aims. In this work, we aim to obtain corresponding constraints in the case of higher dimensional random fields and test them in a 
more realistic context. 

Methods. We develop numerical methods of computing the constraints on correlation functions that are also applicable for two- 
and three-dimensional fields. To test the accuracy of the numerically obtained constraints, we compare them to the analytical results 
for the one-dimensional case. Finally, we compute correlation functions from the halo catalog of the Millennium Simulation, check 
whether they obey the constraints, and examine the performance of the transformation used in the construction of the quasi-Gaussian 
likelihood. 

Results. We find that our numerical methods of computing the constraints are robust and that the correlation functions measured 
from the Millennium Simulation obey them. Even though the measured correlation functions lie well inside the allowed region of 
parameter space, i.e., far away from the boundaries of the allowed volume defined by the constraints, we find strong indications that 
the quasi-Gaussian likelihood yields a substantially more accurate description than the Gaussian one. 

Key words, methods: statistical - cosmological parameters - large-scale structure of the Universe - galaxies: statistics - cosmology: 
miscellaneous 


1. Introduction 


The two-point correlation function £ is been a very common 
tool in cosmology, although an increasing amount of astronom¬ 
ical literature deals with higher order statistics. Whenever cor¬ 
relation function measurements are used in a Bayesian frame¬ 
work to determine cosmological parameters, the probability dis¬ 
tribution function (PDF) of the correlation function is needed. 
Usually, this likelihood, -C(£). is assumed to be a multivariate 
Gaussian distribution (see, for example, an analysis of the corre¬ 
lation function of the cosmic microwave background bv Seliak 
& Bertschinger d 19931) . or co mmo n methods of baryon acoustic 
oscillations detection (e.g., bv lLabatie et alJl2012h l. 


However, the Gaussian approximation of £.(£) is not nec¬ 
essarily well justified in all cases and may not always provide 
the level of precision required from statistical tools that are used 
to analyze state-of-the-art astronomical data, for example, non- 
G aussi anities in the cosmic shear likelihood were detected by 
iHartlan et alJ (12009), In the case of third-order cosmic shear 
statistics, however. Isimon et al.l (1201 5l) recently have found, at 
least in current state-of-the-art surveys, that a Gaussian likeli¬ 
hood is a reasonably good approximation. This agrees with re¬ 
sults for the bispectrum covariance put forward bv lMartin et alJ 
(120121) . As an additional remark, objections against the use of 
Gaussian likelihoods as a ‘safe default’ have been raised in cases 
where one lacks knowledge of the exact form of the likelihood, 
as pointed ou t, f or example, in p ower spectrum analyses by 
(CaiTonl (120131) and lSun et alJ rt2013h . 


A very strong argument against the Gaussianity of Ji(£) 
is the existence of fundamental constraints that stem from the 
non-negativity of the power spectrum and was put forward by 


[Schneider & Hartlapl (l2009i ). hereafter SH2009. That correla¬ 
tion functions cannot take arbitrary values immediately implies 
that the Gaussian approximation cannot be fully correct, since a 
Gaussian distribution has infinite support. To remedy this, one 
might be tempted to use a Gaussian likelihood for £ and include 
the constraints by simply incorporating priors that are zero out¬ 
side the allowed region. However, as shown in previous work 
(see Figs. 4 and 5 in SH2009), the shape of the distributions 
of £ are strongly affected by the constraints, even well inside 
the admissible range and thus a more comprehensive solution is 
needed. 

Of course, it would be preferable to obtain the true PDF of 
£ analytically, which is feasible only for the uni- and bivariate 
cases, even assuming one-dimensional Gaussian random fields, 
as shown bv [Keitel & Schneider ( 201 1). Their results are a cru- 
cial ingredient of the quasi -Gaussian approach introduced in 
IWilking & Schneider! (l2013h . hereafter WS2013. There, we use 
the aforementioned constraints to transform the correlation func¬ 
tion into an unconstrained variable, where the Gaussian approx¬ 
imation is expected to hold to higher accuracy. Using numerical 
simulations, we show that for the correlation functions of one¬ 
dimensional Gaussian fields, this ‘quasi-Gaussian transforma¬ 
tion’ performs very well, meaning that it transforms £ into a vari¬ 
able that is highly G aussian. When we make use of the analytical 
univariate p( f) from iKeitel & Schneider ( 2011 ). this transforma¬ 
tion can then be exploited to construct the quasi-Gaussian likeli¬ 
hood for £. As presented in WS2013, the new description of -£(£) 
agrees well with the distributions obtained from simulations and 
has an impact on the results of Bayesian parameter estimation, 
as shown in their toy-model analysis. 
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To date, a major caveat of the quasi-Gaussian approach stems 
from the fact that the analytical computation of the constraints 
presented in SH2009 is only optimal for one-dimensional ran¬ 
dom fields. This severely limits the set of possible applications of 
the results presented in WS2013. In ISect. 21 we develop numeri¬ 
cal methods to compute the constraints on correlation functions 
that are also applicable to higher dimensional random fields, 
check their robustness, and compare the numerically obtained 
constraints to the analytical results for the one-dimensional case. 
In ISect. 3l we then apply the derived methods in an astrophysi- 
cal context, i.e., to correlation functions measured from the halo 
catalog of the Millennium Simulation. We discuss some practi¬ 
cal aspects of measuring /' and show that the correlation func¬ 
tions obtained from the simulation clearly obey the constraints. 
Furthermore, we examine the performance of the quasi-Gaussian 
transformation: By comparing the skewness and kurtosis of the 
transformed and the untransformed correlation functions, we ar¬ 
gue that the quasi-Gaussian PDF is a better description of the 
likelihood of correlation functions than the Gaussian one. We 
conclude with a brief summary and outlook in ISect. 4l 


2. Numerical computation of the constraints on 
correlation functions 


We consider the two-point correlation function of a random field 
g(x), which is defined as ( : (x.y) = (g(x) g*(y)). It is related to 
the power spectrum via Fourier transform. If assuming isotropy, 
this can be written as 

/ d"k C d k 

~^—^P(\k\) exp(i& ■ x) = J —k n ~ l P(k)Z n (ks), (1) 

where s = |jc|, and the dimensionality n of the underlying random 
field determines the function Z„(rf). For a one-dimensional field, 
becomes a cosine transform; in the 2D case, Z 2 (?/) = 
the Bessel function of the first kind of zero order; and for 
a 3D random field, Z^tj) = jo(rj) is the spherical Bessel function 
of zero order. 

As SH2009 show, correlation functions obey fundamental 
constraints, which arise from the non-negativity of the power 
spectrum and are best expressed in terms of the correlation coef¬ 
ficients r„ = )//(()). As it turns out, the constraints can then 

be written in the form 


Eg. (1) 
Join) is 


r,i\(r\,r 2 ,..., r„_i) < r„ < r„ u (n,r 2 ,..., r„_i), 


( 2 ) 


meaning that the upper and lower boundaries on r n are functions 
of the r, with i < n. 

SH2009 use the fact that the covariance matrix C, ; = 
(gig*) = (where g,- = g(i Ax) for a one-dimensional ran¬ 
dom field evaluated at discrete grid points) has to be positive 
semi-definite, to explicitly calculate the constraints in the case 
of homogeneous, isotropic random fields, and show that the con¬ 
straints they obtain are optimal for a one-dimensional random 
field, meaning that no stricter bounds exist for a general power 
spectrum. For higher dimensional fields, the bounds found for 
the one-dimensional case are still obeyed; however, owing to 
the isotropy of the field and the multidimensional integration in 


Eq. (1) tighter constraints hold that have to be computed numer¬ 
ically. 

The procedure to obtain the optimal constraints numerically 
is outlined in SH2009. Rewriting |Eq. (1) and applying a quadra¬ 
ture formula for the integral yields 


r(s) = £(s)/m = J] Vj Z n (kjs), 


(3) 


j= i 


where the coefficients fulfill 0 < l / ; < I and £ Vj = 1. We note 
that this approximation becomes arbitrarily accurate as K —> oo. 

When measuring correlation coefficients for N different sep¬ 
arations Si, each point r = (i q, r 2 ,..., r^), with r,- = £(s,)/£(0), 
in this /V-dimensional space can be written as a weighted sum 
along the curve c(A) = ( Z n (As \),... ,Z„(Asn)), where we used a 
continuous variable A with 0 < A < oo instead of discrete wave 
numbers kf. 


r = X V J e(Aj). 


(4) 


i=i 


Since 0 < Vj < 1 and £ Vj = 1, each point r has to lie within 
the convex envelope of the curve c(A), which corresponds to 
the constraints on the correlation coefficients; for example, by 
constructing the convex envelope of the curve c(A) for two lags 
(ri, r 2 ) in the one-dimensional case, this reproduces the analyti¬ 
cally known bounds r 2u i(ri). 

As a result, to find the constraints only requires describ¬ 
ing the convex envelope of the curve c(A). Unfortunately, there 
does not seem to be a general analytical solution for this prob¬ 
lem, which means resorting to numerical methods; for ex¬ 
ample, the qhull algorithm dBarber et ali 119961 available at 
http: //www. qhull. org) provides an efficient implementation 
for computing, among other things, convex hulls. It is, however, 
limited to inputs of dimensionality lower than nine, meaning that 
it is only applicable for a maximum number of separations of 
N = 8. Although this is not a requirement for computing the 
constraints, we use equidistant lags throughout this work, denot¬ 
ing s„ = n As. 


As an example of determining the constraints. Fig. 1 shows 
the curve c(A) in the r\ - r 2 -plane, plotted in black up to as 
high as A — 50 for illustrative purposes, as well as its convex 
hull. For a given r \, the upper and lower bounds on r 2 are given 
as intersections with the red hull. This method can, of course, 
be generalized to higher dimensions, e.g. the determination of 
r Su | from r\, . . r 4 , where the convex hull is a hypersurface in 

a five-dimensional space. Following this procedure, we devel¬ 
oped a code to compute the constraints for one-, two-, and three- 
dimensional fields that we use in our analysis of correlation func¬ 
tions measured from the Millennium Simulation in ISect. 3l 


2.1. Comparison of the constraints for one-dimensional fields 

Below, we test the numerical method to obtain the constraints. To 
do so, we compare the numerical results to the analytically com¬ 
puted bounds. Since an analytical calculation of the constraints is 
only possible for one-dimensional random fields and equidistant 
lags, we limit our comparison to this case. Throughout this work, 
we use a gridded approach and denote = %(s n ) = %(n As), 
where As = L/N is the separation between adjacent grid points, 
and L denotes the length of the field. 

There are several ways of testing our methods of comput¬ 
ing the constraints: Most straightforward is to compare the con¬ 
straints from the two methods directly, i.e., to compute the up¬ 
per and lower bounds r nu and r n \ both analytically and numeri¬ 
cally, and check how much they differ. An alternative approach 
involves the quasi-Gaussian transformation r„ —> y„, which, as 
briefly explained in ISect. ll is a central ingredient of the quasi- 
Gaussian approximation for the likelihood of the correlation 
functions introduced in WS2013: 

, 2 r n r nu r n \ 

y n = atanh-. (5) 

F mi Yn\ 
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Fig. 1. Example of the curve c(A ) for a two-dimensional random 
field in the r\ - ri-plane, where c(/l) = (Jo(A), Jq(2A)). The red 
circles and the line connecting them show the convex hulls de¬ 
termined by qhull. 


Since this transformation is the main application for the con¬ 
straints, it can - and should be - applied as a means to com¬ 
pare the analytically and numerically obtained bounds, namely 
by using the different sets of constraints in the transformation 
and comparing the resulting y n . 

As previously described, the constraints on r„ are functions 
of the correlation coefficients with lower lags, and as such, we 
need input values for r \,..., i to be able to compute and com¬ 
pare the different r n \ and r na . Again, two possibilities arise: To 
provide input values that are close to ‘real-life’ applications, we 
can use realizations of correlation coefficients obtained from nu¬ 
merical simulations (see WS2013 for an efficient way to gener¬ 
ate realizations of the correlation function of a one-dimensional 
Gaussian random field). However, this obviously requires as¬ 
sumptions about the underlying random field and, in particular, 
its power spectrum. Consequently, a more general approach is 
to draw the input correlation coefficients for computing the con¬ 
straints randomly, i.e., from a uniform distribution over the al¬ 
lowed range, r n e ]r n \, r na [. Due to the nature of the constraints, 
this is an iterative procedure, meaning that one has to draw 
r\ € ]rn, n u [, compute r 2 i ?u from this ri, then draw 7-2 e ]r 2 i, r 2 U [ 
to determine r^u, and so on. 

A comparison of the analytically and numerically obtained 
bounds r„ u and r n \ is shown in Fig. 2 where we plot the dif¬ 
ferences /- ar ' a - r" 1 " 11 as functions of n. For each bound r, m 1 , the 
required input values of the correlation coefficients r, with i < n 
are drawn uniformly, as previously described. To perform a sta¬ 
tistically significant check, this procedure is repeated 500 times, 
meaning that we generate 500 realizations of the input correla¬ 
tion coefficients and compute the upper and lower bounds both 
numerically and analytically for each realization. The values 
plotted in the figure are obtained by averaging the difference 
between the analytical and the numerical values over the 500 
realizations. In addition, we also investigate how much impact 
the sampling of the convex hull of the curve c(A) has on the ac¬ 
curacy of the numerical bounds. It turns out that it is sufficient in 


Fig. 2. Difference between the analytically and numerically ob¬ 
tained bounds, averaged over 500 realizations. The upper three 
sets of points correspond to the difference r ana - r" um , whereas 
the ones with negative values show r ana - r™" 1 . Furthermore, the 
different symbols denote the number of steps used to sample the 
convex hull of the curve c(/l) for values of 0 < A < 2n, namely 
100 (red crosses), 200 (blue circles), and 300 (green triangles) 
steps. 


all cases to sample the curve c(A) for values of 0 < A < 2n, since 
going to higher values of A has no impact on the volume within 
the convex hull because of the periodicity of the function Z n (ks) 
in Eq. (1) By way of comparison, we vary the sampling rate of 


the convex hull; i.e., we sample it using 100 (red crosses), 200 
(blue circles), and 300 (green triangles) steps. The upper three 
sets of points show the difference r ana - r™ m between the analyt¬ 
ical and numerical upper bounds, whereas the lower ones depict 


Three conclusions can be drawn from Fig. 2 First, the de¬ 
viation of the numerically obtained bounds from the analytical 
ones shows a tendency to grow with n, which is to be expected, 
since the sampling of the convex hull becomes more challenging 
with increasing dimensionality. (The numerical and analytical 
bounds on r\ do not differ at all, since ri u j = ±1.) Second, the 
impact of this sampling has a strong impact on the accuracy of 
the numerical calculation of the bounds; namely, the difference 
between the numerical and the analytical results decreases by 
about a factor of three when doubling the number of steps used 
for the convex hull sampling. Actually, this sampling is the lim¬ 
iting factor for the accuracy of the numerical bounds, as can be 
seen from the third observation: In the case of the upper bounds, 
the numerical results are systematically smaller than the exact 
analytical values, whereas for the lower bounds, the numerical 
values are too high. This effect is an expected consequence of 
the non-continuous approximation for the smooth hull; as a re¬ 
sult of convexity, the hyperplanes that link the points used to 
sample the hull always have to be located inside the hull. 

In summary, the accuracy of the numerical constraints can be 
increased by improving the sampling of the hull. While a larger 
number of steps would, presumably, improve the results even 
further, using more than 300 steps for the sampling becomes im- 
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Fig. 3. Difference between the y n computed using the analyti¬ 
cally and the numerically obtained bounds, averaged over 500 
realizations, with error bars showing the standard deviations. 
The input values for the computation of the bounds used in the 
transformation r„ —> y„ stem from 500 simulated realizations of 
the correlation function on a one-dimensional Gaussian field of 
length L with N — 32 grid points and a Gaussian power spectrum 
of width ko , with Lk {) = 80. In the case of the black, solid points, 
the convex hull of the curve c(A) is sampled using 100 points 
for the range 0 < A < 2n. By way of comparison, we show the 
corresponding results for sampling rates of 200 (red cross) and 
300 (blue circle) in the case n = 8. 


practical because of the computational costs. However, as the 
following tests demonstrate, using 300 steps seems sufficiently 
accurate. 

As mentioned before, another important check for the ac¬ 
curacy of the numerical methods that are used to compute the 
bounds is to apply the quasi-Gaussian transformation r n —> y„ 
and to compare the resulting y n . Below, we adopt correlation 
coefficients from simulations instead of uniformly drawn ones 
as input for the computation of the bounds. In particular, we use 
500 realizations of the correlation function on a one-dimensional 
Gaussian field of length L with N = 32 grid points and a 
Gaussian power spectrum, where the field length and the width 
of the power spectrum are related by Lko = 80. (For details on 
the simulations used in this work, we refer to WS2013.) 


For each simulated realization of the correlation coefficients, 
we compute the bounds r„ u and r n \ for each n, both numerically 
and analytically, and use them to transform r n to y„, as defined 



the error bars denote the standard deviations. For the sake of 
clarity, we only show the impact of the number of steps used in 
the convex hull sampling for the case n = 8, where the standard 
deviations are largest. 

The accuracy of the numerical approximation again shows 
a strong dependence on the number of steps used to sample the 
convex hull. Nevertheless, the difference between the values of 
y n , computed using the analytical and the numerical bounds, be¬ 


comes very small when using 300 steps. As a result, we con¬ 
clude that the problem of the non-continuous approximation of 
the curve c(A) and its convex hull can be tackled and that using 
300 steps in the sampling yields sufficiently accurate bounds. 


3. Application to the Millennium Simulation 

So far, our studies about the constraints on correlation func¬ 
tions and the quasi-Gaussian likelihood have been performed 
in a general, mathematical framework. In this section, we in¬ 
vestigate our results in a more astrophysical context by apply¬ 
ing them to cosmological correlation functions measured in the 
Millennium Simulation (ISpringel et alj[2005h . Thus, we aim to 
check the relevance of the constraints that originally stem from 
purely mathematical properties, since they are based on the fact 
that £ is the Fourier transform of a positive quantity (the power 
spectrum), in a more practical context, where £ is measured us¬ 
ing an estimator. The size of the Millennium Simulation enables 
us to easily measure multiple realizations of £, thereby provid¬ 
ing an approximate determination of the underlying probability 
distribution and, consequently, a statistical analysis. 


3.1. Computing correlation functions 


Below, we compute the correlation function of dark matter halos 
in the Millennium Simulation. Because we are not interested in 
redshift evolution, we only use the halo catalog from the z - 0 
simulation snapshot, from which we then select typical galaxy- 
mass halos by choosing a mass cut > 10 12 h~ l M e , which 
yields a total number of ~ 440 000 halos. However, to per¬ 
form a statistical analysis, we require different realizations of 
the correlation function. For this reason, as a first attempt, we 
divide the full simulation cube into 1000 subcubes of volume 
50 3 [hr 1 Mpc) and measure £ in each of the subcubes. 

Along with the halo catalog from the simulation we also 
need a random catalog, so, for each subcube, we draw halo posi¬ 
tions uniformly. We then determine the number of halo pairs for 
given pair separations in both the data and the random catalog, as 
well as the cross-correlation. From the count rates DD(s), RR(s), 
DR(s) (normalized to account for different numbers of halos in 
both the random catalog and the halo catalog from each subcube) 
at different pair separations s , we compute £ using an estima¬ 
tor. While Landy-Szalay (LS) is the most widely used estimator 
dLandv & Szalav i ; 1993 ). we also aim to test the impact of the 
choice of estimator on the con straints, and adapt the follo wing 
common set of estimators from lVargas-Magana et al.l (1 2013 ): 

£ph = - 1. Peebles_&Hauser Q974); 

= T Hewett (1982) : 

£ DP = ££L-l, Davis & Peebles (1983) : 


£ H = 22|P-1, Hamilton (1993) : 

£ ls = dd-idr+rr t Landv & Szalav (1993) . 

As mentioned in ISect. 21 to calculate and test the constraints, 
we measure the correlation function at equidistant lags, i.e., 
£„ = £(n ■ As), where the maximum number of lags is n = 8, 
an effect of the limitations of the numerical computation of the 
constraints. 

The size of the random catalog also merits some discus¬ 
sion: Ideally, it should be infinitely large, i.e., N mn d —> °°- 
However, the computation of the pair separations is the most 
time-consuming step in the calculation of £ and, consequently, 
the number of halos in the random catalogs for each of the 
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Fig. 4. Correlation function from 1000 subcubes of the 
Millennium Simulation, computed using the LS estimator. The 
points and error bars show the correlation function • As) 

for As = 5 li 1 Mpc (see text for details) averaged over the 1000 
subcubes of side length 50 Ir 1 Mpc, as well as the standard de¬ 
viation. For the blue circles, the random catalog for each sub¬ 
cube contains 100 halos, as opposed to 10000 halos for the red 
crosses. 


1000 subcubes is subject to practical limitations. We study the 
impact of the random catalog size in Fig. 4 Here, we show 


the correlation function for an exemplary choice of lags with 
As = 5 h~ l Mpc, i.e., we measure ■■•,^8 at l a g s of 

5,10,.. .40 h 1 Mpc. In practice, we need to allow for a range of 
pair separations to obtain sufficiently large numbers of pairs. To 
do this, we adapt a bin size of width 1 h~ l Mpc, so, for example, 
in the computation of £\ we use all pairs with separations ranging 
from 4.5 to 5.5 li 1 Mpc. For the auto-correlation fo, i.e., the cor¬ 
relation function at zero lag (which we do not plot in the figure, 
but which is required for the calculation of the constraints), we 
count all pairs with very small separations, e.g., s < 1 hr 1 Mpc. 
(We refer to the next section for a discussion on the measure¬ 
ment of A) and on the choice of lags.) The points and error bars 
show the mean and standard deviation over the 1000 subcubes 
of side length 50 h 1 Mpc; we use the LS estimator, which has 
been shown to be less sensitive to the size of the random cata¬ 
log than others (see lKerscher et al .1-2000) . For the blue circles, a 
small random catalog (Wand = 100 halos for each subcube) was 
used, whereas the choice of Wand = 10000 for the red crosses 
results in noticeably smaller standard deviations over the 1000 
realizations, at the cost of a longer computation time. 

Hence, we aim to find a trade-off between those two val¬ 
ues. First, although the mean of the correlation functions for the 
two random catalog sizes plotted in Fig. 4 do not seem to dif¬ 
fer very much at first sight, choosing the catalog size as small 
as Wand = 100 is a quite extreme case. This is because a large 
fraction of the realizations yield a diverging auto-correlation <fo 
as a result of the count rates in RR or DR being zero, at least 
when measuring £q as previously described. However, when in¬ 
creasing the random catalog to 1000 halos per subcube, the mean 
correlation function for non-zero lag shows a deviation of only 


Fig. 5. Histograms of the distributions of overdensities in the 
subcubes of the simulation volume, where e denotes the num¬ 
ber overdensity of halos in the subcube relative to the mean halo 
number density in the whole simulation box, see Eq. (6) The 
colors indicate the number of subcubes the simulation box was 
divided into. 


~ 1 % compared to the mean (■' for Wand = 10000 (and even 
here, about a tenth of the realizations show a diverging <;<)). At 
Wand = 5000, no such problems occur, while even the error bars, 
as plotted in the figure, become indistinguishable from those at 
Wand = 10000. This means that a random catalog size of 5000 is 
a reasonable trade-off between accuracy and computational ex¬ 
penses. 


An additional observation from Fig. 4 is that £ becomes neg¬ 
ative for larger lags, i.e., around 20 - 25 h -i Mpc. The reason for 
this is an integral constraint (see, for example, Landv & Szalav 
1993 1 that arises when measuring correlation functions in finite 
volumes, where the global mean density is unknown and is usu¬ 
ally approximated by the mean observed density. 

In our case, one way to assess this issue is to decrease the 
number of subcubes in our analysis, and make them larger and 
more representative for the whole box, while at the same time 
measuring £ at the same lags as before. Beside lessening the im¬ 
pact of the integral constraint on small lags, this has two addi¬ 
tional effects: First, with fewer subcubes, we obtain fewer real¬ 
izations of which can pose a challenge for a statistical analysis, 
and second, as the number of halo pairs per cube becomes larger, 
the scatter over the measured realizations of £ decreases. To de¬ 
cide on the number of subcubes necessary to make the subcubes 
as representative as possible for the whole simulation volume, 
we estimate the overdensity e in each subcube by comparing the 
mean number density in the subcube to the one from the whole 
simulation volume, using 


^sub — ^box (1 "k €”) 5 


( 6 ) 


and examine the distributions of e by plotting them as histograms 
Here, we slice the simulation volume into different 


in 


Fig. 5 


numbers of subcubes and compute the overdensity in each sub¬ 
cube. It is clear that the distribution p(e) is very broad for the 
value Wubcubes = 10 3 used so far, and, as expected, it becomes 


5 





































P. Wilking et al.: Constrained correlation functions from the Millennium Simulation 


m 


S3 

fJLP 


in 

o 


o 

o 


N su b cu b es — 8 
I N su bcubes — ^ 


~r~ 

4 


-r - 

7 


n 

Fig. 6. Correlation function from the Millennium Simulation, 
computed using the LS estimator, a random catalog size of 5000, 
and a lag separation of As — 5 h 1 Mpc. The points and error 
bars show the mean and standard deviation computed over the 
subcubes of the simulation, where the simulation box was sliced 
into 8 1 subcubes for the blue data points, as opposed to 5 3 for 
the red ones. 


quite narrow for the case of 5 3 subcubes, which indicates that 
the integral constraint does not pose a large problem in this case. 
The resulting correlation functions for the cases of 8 3 and 5 3 sub¬ 
cubes (and otherwise the same parameters as before, i.e., same 
lags and random catalog size) are shown in Fig. 6 As can be 
seen, slicing the simulation volume into 5 3 subcubes yields rea¬ 
sonable results, i.e., a non-negative correlation functions with a 
sufficiently small variance. 

Finally, we briefly evaluate the choice of estimator. To do so, 
instead of measuring £ n at a few different lags n, it is advisable to 
compute 4 ; (.v) for all lags in each subcube. Since this is obviously 
not practicable, we compute it at a high number of lags, meaning 
that we divide the range of pair separations into adjacent bins of 
width 0.2 hr 1 Mpc. 

The correlation functions (average over the 5 3 subcubes and 
using a fixed size of random catalog for each subcube, namely 
iV ra nd = 5000) are shown in Fig. 7 Here, the different colored 
lines denote the five estimators, and the gray-shaded region de¬ 
picts the standard deviation over the 125 realizations in the case 
of the most commonly used LS estimator. For clarity, in the left 
panel, we plot scales from 8 to 40 h~ l Mpc, whereas the right 
panel shows the correlation function for very small lags, i.e., 
4-8 lr 1 Mpc. Clearly, the numerous estimators yield very sim¬ 
ilar results, in particular compared to the standard deviation of 
the 125 realizations. 


3.2. Testing the constraints 

In this section, we investigate whether the correlation functions 
computed from the halo catalog of the Millennium Simulation 
obey the numerically obtained constraints. While we make the 
case in ISect. 2TT1 that using 300 points to sample the hull yields 
sufficiently good agreement between the numerical and the ex¬ 


act analytical values of the constraints, we have to restrict our¬ 
selves to 270 points in the case of a 3D random field owing to 
the computational costs. Although the convex hull only has to be 
computed once and can then be used to determine the constraints 
for all sets of correlation coefficients, sampling the hull for a 3D 
random field with the given accuracy poses memory problems 
for the qhull algorithm, and is beyond the scope of this research. 
However, this does not pose a problem: When we compare the 
accuracy of the numerical constraints, as plotted in Fig. 2 it is 


apparent that the improvement in accuracy, when going from 200 
to 300 steps, is far smaller than the one from 100 to 200, and thus 
we expect 270 points to be accurate enough. 

To test the constraints, for each realization we compute the 
correlation coefficients r n = T n /To as well as the upper and lower 
bounds, r nu and r n \. It turns out that the width of the Aj-bin has 
a strong influence, in particular on the width of the distributions 
of the correlation functions r„. For example, we first choose a 
relatively broad bin, i.e., we measure by averaging over all 
pair separations from 0 to 2 li~ l Mpc. This choice is primarily 
motivated by the fact that increasing the spread of the correlation 
coefficients over the 125 realizations allows us to test how close 
to the edges of the allowed region the r n move. Toward the end 
of this section, we study the impact of the width of the £o-bin in 
more detail. 

One question that arises is how to visualize the constraints. 
The simplest approach to this is to use scatter-plots, with dots 
for the individual realizations. An example in the r i - r 2 -plane 
is shown in |Fig. 8| Here, the red dots show the different real¬ 
izations of r\ and r 2 , computed for the subcubes using the LS 
estimator. Additionally, we plot iso-density contours that con¬ 
tain 68, 95, and 99.7 % of these realizations. For the lefthand 
panel, we sliced the simulation volume into 1000 subcubes, as 
opposed to 125 for the righthand panel. As explained in the pre¬ 
vious section, the higher number of subcubes greatly increases 
the spread of the correlation functions, which can also be clearly 
observed in /--space. (Even for the high number of subcubes, the 
integral constraint is expected to be negligible for the correlation 
functions at small lags.) In both panels, the upper and lower blue 
lines represent the constraints, i.e., r 2 u .i(/'i), which we compute 
numerically for each realization of r\ shown in the figure and 
plot as connected lines. All realizations clearly lie well inside 
the constraints, particularly when compared to results for purely 
Gaussian (one-dimensional) random fields with fiducial power 
spectra (see, for example, similar figures in related works, such 
as Fig. 3 from WS2013). 

As an additional way of depicting the constraints, we ap¬ 
ply a part of the quasi-Gaussian transformation from |Eq. (5)| to 
map the allowed range of the correlation coefficients to (-1, +1), 
namely by transforming the correlation coefficients r„ to 


2 r n /'/m r n \ 

x„ = ---. (7) 

bm r n 1 

To visualize the constraints more clearly, we use a modified 
version of box-and-whisker plots, meaning that we display our 
samples {r,,} and {x„} as boxes, where the upper and lower bor¬ 
ders show the first and third quartiles of the sample, i.e., the 
values that split off the upper and lower 25 % of the data. As 
is common practice, we also show the sample median (i.e., the 
second quartile) as a line inside the box, as well as two whiskers. 
Usually the ticks at the end of the whiskers denote the minimum 
and maximum of the data (in the most widely used type of box- 
and-whisker plot). Here, we use them to denote the upper and 
lower constraints: Since r„ Ui 1 are functions of all r, with i < n. 
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s [h -1 Mpc] 


Fig. 7. Correlation function measured for ‘all’ lags (see text for details) as a function of the pair separation s , with a random catalog 
size of iV ra nd = 5000. The lines denote the mean £(s) measured in the 125 subcubes using the different estimators listed in IScct. 3~i~l 
and the gray shaded region shows the standard deviation. In the left panel, the ,v-range from 8 to 40 li~ l Mpc is plotted, and the right 
panel shows the results for very small lags, from 4 to 8 hr 1 Mpc. 




r i n 

Fig. 8. Correlation coefficients r\ and ri measured from the halo catalogs in the subcubes of the Millennium Simulation, using the 
LS estimator, where we slice the simulation volume into 1000 subcubes for the left panel and 125 for the right one, and the random 
catalog for each subcube contains 5000 halos. In both cases, we measure £ at lags of separation A,v = 5 h~ l Mpc, and use all halo 
pairs with pair separations of 0 to 2 h 1 Mpc to compute The red dots show the 1000 (125) realizations, while the black lines 
are iso-density contours that contain the given percentages of the realizations. The upper and lower constraints riyji/'i), computed 
individually for each realization of r \, are shown as blue lines. 


we show the mean r„ u and r n \ over all realizations for plots in r- 
space. For the transformed values x ni the bounds are simply ±1, 
so there is no need to average over the realizations. 


Fig. 9 shows box-and-whisker plots of r„ and x„ at all eight 


lags n, where we use the same lags and random catalog size as 
before, as well as the LS estimator. It can be seen that the con¬ 
straints are clearly obeyed, and although the distributions be¬ 
coming broader for increasing lags, the boxes showing the upper 
and lower quartiles only occupy a small portion of the allowed 


region. The distributions are not necessarily centered within the 
allowed region, which is not surprising, since their exact shape 
and position also depend on the underlying power spectrum. The 
choice of estimator has barely any impact on the widths and po¬ 
sitions of the distributions within the allowed region, which is 
to be expected, since Fig. 7 already illustrates that the different 
estimators yield quite similar results. 


As mentioned at the beginning of this section, the main in¬ 
fluence on the variance of the distributions in A and correspond- 
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n n 

Fig. 9. Box-and-whisker plots for the correlation coefficients r n and the transformed values x„, as defined in |Eq. (7)] where each data 
point shows the upper and lower quartile (edges of the box), median (line inside the box), and mean upper and lower boundaries 
(whiskers) for the 125 realizations measured from the simulation subcubes. These use the same lags, random catalog size, and 
estimator as before (see text and previous figure captions for details). 




As; bin width; ^o - bin width [h 1 Mpc] As; bin width; ^o - bin width [h 1 Mpc] 

Fig. 10. Box-and-whisker plots of the transformed correlation coefficients at smallest and largest lag for varying lag separation and 
bin width. The triple labeling of each distribution gives the lag separation As, the bin width of pair separations at non-zero lag (i.e., 
for £i... £s), and the width of the 4 r o-bin. For example, for the second case shown in each panel, we measure <fo from halo pairs with 
separations from 0 to 1 h~ 1 Mpc, £\ from those with separations from 4 to 6 lr ] Mpc, £2 from 9 to 11 h 1 Mpc, and so on. 


ingly in r- and x-space, seems to be the width of the t'o-bin. In 
|Fig. 10| we investigate this observation and also study the im¬ 
pact of the choice of the separation between the lags at which 
we measure £. In the two panels of the figure, we show box- 
and-whisker plots of the transformed correlation coefficients at 
smallest and largest lag, i.e., xi and xs, and we vary the sepa¬ 
ration As of the lags as well as the bin widths of the pair sepa¬ 
rations used to measure £0 and the correlation functions at non¬ 
zero lag, £1 . ■ - £'h- In the case of the four left-most distributions 
in each panel, we use a lag separation of As = 5 h~ l Mpc, 


where we adapt a bin width of 1 h 1 Mpc for £\ ... t'x for the 
first and second distribution, and a bin width of 2 /T 1 Mpc for 
the third and fourth ones. In both cases, we separately use a nar¬ 
row and a broad bin width for the measurement of fo (also 1 and 
2 h 1 Mpc). The figure illustrates that the width of the distribu¬ 
tions of x„ is mainly determined by the Arbin size, whereas the 
width of the bins for £ n at lags n > 0 barely has any influence. 
This is due to the structure of the quasi-Gaussian transforma¬ 
tion, where 4 ; o appears in every correlation coefficient r „, and, 
as a result, in the computation of every lower and upper bound. 
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Fig. 11. Test for the univariate Gaussianity of the {£}-and {y}-samples obtained from the Millennium Simulation, using a lag separa¬ 
tion of As = 5 h 1 Mpc and bin widths of 1 h 1 Mpc for all including £o. The black circles and red triangles show the univariate 
skewness and kurtosis of the distributions p(£ n ) and p(y n ), computed over 125 (solid curves) and 1000 (dashed curves) subcubes of 
the simulation volume. For the blue curves, we draw 100 Gaussian samples with the same mean, covariance matrix and sample size 
as the distributions p(y„) in the case of 125 subcubes, and plot the mean and standard deviation of their skewness and kurtosis. 


The impact of the width of the £o-bin is particularly strong for 
small-lag distributions, which it also shifts, as can be seen from 
the distribution of x\. In particular, this shift is larger than a case 
where we measure at different lags altogether, as illustrated 
for a lag separation of A,y = 3 h 1 Mpc in the fifth distribu¬ 
tion shown in the figure. In this context, it is important to stress 
that the problem of how to measure 4 ; (l in practice is well-known 
since, in most applications, it is difficult to measure £ at very 
small lags. As we show, however, this poses a particularly diffi¬ 
cult challenge when analyzing measured correlation functions in 
a quasi-Gaussian framework, since here, the exact determination 
of fo is vital; the auto-correlation function enters everywhere, 
because one would always transform f to y (or at least to r) for 
an analysis involving the constraints. 

In summary, all correlation functions we measured from the 
Millennium Simulation are quite far away from the edge of the 
allowed region. This finding seems to hold irrespective of the 
choice of estimator, lags, etc., providing that £ is measured in a 
‘sensible’ way. As an example, using very small random catalog 
sizes does indeed yield single realizations outside the allowed 
region. 


3.3. Quality of the Gaussian approximation in g andy-space 


In this section, we use the correlation function samples measured 
from the Millennium Simulation to assess the quality of a quasi- 
Gaussian approach. Similar to the tests shown in WS2013 for 
simulated correlation functions, we transform £ to y as defined 
in 


£ 


Eq. (5) and test the Gaussianity of the distributions in y and 


recause the Gaussianity in y-space is a central prerequisite for 


the accuracy of the quasi-Gaussian likelihood. 


While it would be preferable to assess the quality of the 
quasi-Gaussian approximation directly, i.e., to check how well 
the quasi-Gaussian PDF agrees with p(£), as obtained from the 
Millennium Simulation, computing the quasi-Gaussian PDF still 


requires measuring the underlying power spectrum, which is be¬ 
yond the scope of this work. In real life, however, one would 
usually transform the measured correlation function to y-space 
to perform a Bayesian analysis and, thus, the Gaussianity of 
p(y ) is pivotal. Even so, knowledge about the underlying power 
spectrum would still be required to make use of the analytically 
known p(£;o). 

In the literature, various tests for Gaussianity exist. In this 
study, we focus on the calculation of moments; in particular, we 
compute the skewness and kurtosis, which are defined in such a 
way that they are zero for a Gaussian distribution. In the univari¬ 
ate case, the skewness y of a distribution p(x) reads 


(x - p) 3 \ _ m 3 

rr 3 I - 3/2 ’ 

<r I m 2 


( 8 ) 


where m ,■ = {(x-p)') denotes the central ;th-order moment. Thus, 
y is essentially the (renormalized) third-order moment, and the 
kurtosis 


k — 


(x -p) 4 \ m A 

a 4 I ml 


3 


(9) 


is closely related to the fourth-order moment. In the multivariate 
case, we use the definitions established bv lMardia (1970, 1924), 
who define the skewness of a d-v ariate distribution as 


** = b Z Z { ( * ; - ^ c ~ 1 ( Xj ~ ^)) 3 ’ 

/=! J= 1 


( 10 ) 


where n is the sample size, and p and C are the sample mean and 
covariance matrix. The kurtosis measure reads 


| n ^ 

Kd = - 'Yj {(*'■ - P) T C~ l (Xj - p)}~ -d(d + 2), (11) 

1=1 

where we subtract the last term to ensure that a perfectly 
Gaussian sample yields Kd = 0. 
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Fig. 12. Multivariate Mardia’s skewness and kurtosis of the n-variate distributions of the (f')-and {y}-samples, obtained from the 
Millennium Simulation, and of corresponding Gaussian samples, using the same parameters as before (see previous figure caption). 
In contrast to the previous figure, we plot two curves for Gaussian samples: For the solid (dashed) curve, we draw Gaussian samples 
with the same mean, covariance matrix, and sample size as the corresponding distributions in y-space for the case of 125 (1000) 
subcubes; the blue squares and error bars show the mean and standard deviation of the skewness and kurtosis of the 100 samples. 


To test the impact of the quasi-Gaussian transformation on 
Gaussianity, we transform each realization of the correlation 
function (measured for eight lags of separation A,v = 5 h 1 Mpc 
with bins of width 1 h 1 Mpc for £o ... £h) to y and compute 
skewness and kurtosis of the distributions in and y-space. 
Analogously to our tests in WS2013, we also draw Gaussian 
samples with the same mean and covariance matrix as our sam¬ 
ples jy), both for comparison and to account for small sample 
sizes. The results for the univariate distributions are plotted in 
Here, we show the skewness and kurtosis of the distri- 
p{% o). ■ ■ •, P(ft) and p(yi), p(y g ); for the solid lines, 
we sliced the simulation volume into 125 subcubes, whereas the 
dashed lines correspond to the case of 1000 subcubes. By way 
of comparison, the blue curves show the skewness and kurtosis 
of corresponding Gaussian samples, where for the sake of clar¬ 
ity, we only plot the curves for 125 subcubes. Because skewness 
and kurtosis fluctuate quite significantly for this small sample 
size, we draw 100 Gaussian samples of size 125 and compute the 
skewness and kurtosis of each sample; the blue squares and error 
bars show the mean and standard deviation of the skewness and 
kurtosis of the 100 samples. We only include the purely Gaussian 
samples for comparison and to check how close to zero the mea¬ 
sured skewness and kurtosis are for these small sample sizes. 
Specifically, they are not meant to provide any insight into the 
absolute scale of the non-Gaussianity observed in £ : . Evidently, 
the distributions in y are far more Gaussian than those in £ (with 
the exception of p(ifo) in the case of 125 subcubes) and, of partic¬ 
ular note, show a kurtosis comparable to the Gaussian samples. 


Fig. 11 
butions 


Since the Gaussianity of the univariate distributions does 
not imply Gaussianity of the multivariate PDFs, we also com¬ 
pute the moments of the n-variate distributions p(go ,... ,£„-i), 
p(yu ..., y„) and of corresponding multivariate Gaussian sam¬ 
ples, and plot them as functions of n , as in Fig. 12 Here, we 


show the results for corresponding Gaussian samples for both 
125 and 1000 subcubes, where the plotted values and error bars 


are the mean and standard deviation of the skewness and kurto¬ 
sis computed over 100 Gaussian samples. While the multivariate 
moments of the Gaussian samples of size 125 are not consistent 
with zero, this is indeed the case for the larger sample size of 
1000. For the dashed curves, i.e., the case of 1000 subcubes, in 
the case of higher n, the integral constraint has a non-negligible 
impact on the measured correlation functions, as explained in the 
previous section. As a consequence, the corresponding skewness 
and kurtosis results should only be considered quantitatively. 
Nevertheless, it is clear that the difference between the level of 
Gaussianity in £- and y-space becomes even larger for the mul¬ 
tivariate case, reaching about one order of magnitude in y and k 
in the case of 125 subcubes. 


As we demonstrated in the previous section, the width of 
the fo-bin, i.e., the range of pair separations used to measure 
the auto-correlation function, has an impact on the distribu¬ 
tions of the correlation coefficients, and thus on those of the 
y„. Accordingly, we vary the £'o-bin width and again study the 
multivariate moments of the corresponding distributions. |Fig. 13] 
shows a similar plot to Fig. 12 but we use a fo-bin width of 2 in¬ 
stead of 1 h 1 Mpc. As it turns out, this yields distributions in y- 
space which are almost perfectly Gaussian, at least in the case of 
125 subcubes, where their moments are hardly distinguishable 
from those of the corresponding Gaussian samples with same 
sample size. As before, it seems that the width of the £o-bin has 
a far higher impact on the results than the bin widths for £) ... £g. 
Actually, using bins of 2 h 1 Mpc for the higher-lag correlation 
functions barely influences the outcome. 


In summary, all tests shown in this section indicate that the 
distributions in y-space are far more Gaussian than those in 
and in some cases even have skewness and kurtosis comparable 
to those of Gaussian samples of the same size. This demonstrates 
the validity of the quasi-Gaussian approach independent of the 
specific parameters used to measure the correlation function. 
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Fig. 13. Multivariate skewness and kurtosis of the (f')-and (y {-samples obtained from the Millennium Simulation and of correspond¬ 
ing Gaussian samples. When compared to the previous figure, we adapt a broader fo-bin, he., we measure the auto-correlation 
function from all halo pairs with pair separations from 0 to 2 h 1 Mpc. 


4. Conclusions and outlook 

Building on SH2009, we have developed numerical methods to 
compute the fundamental constraints on correlation functions. 
We have shown these methods, which are applicable also in the 
case of two- and three-dimensional random fields, to be robust 
and precise, since the numerical computation of the constraints 
for the one-dimensional case reproduces the analytically known 
bounds. We then applied our results to samples of correlation 
functions measured from the halo catalog of the Millennium 
Simulation. After discussing some challenges in the measure¬ 
ment of f, such as the choice of random catalog size and lag 
separation, as well as the question of how to overcome the in¬ 
tegral constraint, we have shown that the correlation functions 
measured from the simulation very clearly obey the constraints. 
Even though all measured correlation functions lie far away 
from the edges of the allowed region, we have demonstrated that 
the quasi-Gaussian quantity y yields significantly smaller non- 
Gaussian signatures (i.e., skewness and kurtosis) than the orig¬ 
inal correlation function £, giving further support to the claim 
that the quasi-Gaussian approximation for the correlation func¬ 
tion likelihood, introduced in WS2013, is a far better description 
than the Gaussian one. 

As a brief outlook on possible future work, one vital im¬ 
provements would be to bypass the current limitation to eight 
lags in the numerical computation of the constraints, since mod¬ 
ern astronomical observations usually measure £ at far more 
lags. Furthermore, the performance of the quasi-Gaussian likeli¬ 
hood in the three-dimensional case should be assessed and com¬ 
pared to the classical Gaussian approach. While this would be 
testable on the samples of correlation functions measured from 
the Millennium Simulation, the most significant advance would 
be the application of our methods to real data, and an investiga¬ 
tion of their impact on cosmological parameter estimation. Aside 
from the current limitation to only eight lags, this would pose 
additional challenges, depending on the area of application: In 
the case of a redshift survey, for example, different constraints 
on the correlation function, measured along and perpendicular 


to the line-of-sight, would hold as a result of redshift space dis¬ 
tortions. Nonetheless, the constraints on correlation functions of 
three-dimensional random fields are, in principle, treatable, de¬ 
spite open challenges and room for improvements. Taking these 
things into consideration, this work opens up a vast field of ap¬ 
plications where Gaussian likelihoods for £ have previously been 
used. 
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